%% Load data
table1 = readtable("I:\MOSCHINI\Lexus\Smith\Crop Rotation\Crop Rotation\CSB\csbDataLong.csv");
table1 = table1(table1.t>=9,:);

%% Breakdown table by LRR and field size and save as .mat files
% Set up variables needed
LRRs = ["F" "G" "H" "K" "L" "M" "N"];
sizegroups = "small";

%% Load estimation table
estimates = readtable("I:\MOSCHINI\Lexus\Smith\Crop Rotation\Crop Rotation\CSB\Estimation results_final.xlsx","Sheet","mlogit parameters");
i = 0;
holder = struct;
for j=1:length(LRRs)
    lrr = LRRs(j);
    for k=1:length(sizegroups)
        size = sizegroups(k);
        i = i+1;
        holder(i).lrr = lrr;
        holder(i).size = size;
        holder(i).vce = readmatrix("I:\MOSCHINI\Lexus\Smith\Crop Rotation\Crop Rotation\CSB\Model\Estimates\vce_final_"+lrr+"small.xlsx","Sheet","Sheet1","Range","C3:NB366");
    end
end


%% Do stuff
%--- Wait bar to show status ---%
tic
f = waitbar(0, "Starting first iteration");

i = 0;
N = 100;
for j=1:length(LRRs)
    lrr = LRRs(j);
    
    for k=1:length(sizegroups)
        %%
        i = i+1;
        size = sizegroups(k);
        
        %--- Get paramter values for the LRR and field size group
        % Logic to get estimates for the group
        group = strcmp(estimates.LRR,lrr) & strcmp(estimates.fieldgroup,size);
        % Get table of estimates from this group
        est = estimates(group,:);
        % Get numerical arrays of b (parameter estimates)
        b = est{:,"b"};
        V = holder(i).vce;
        V = V(1:length(b),1:length(b));
        
        % Filter table to just the current LRR and size group
        %--- Get rows of the data table for the LRR and field size group
        if size=="large"
            rows = strcmp(table1.lrrsym,lrr) & table1.acres>=10 & table1.sameCrop==0;
        else
            rows = strcmp(table1.lrrsym,lrr) & table1.acres<10 & table1.sameCrop==0;
        end
        if strcmp(lrr,"G")
            % Remove regions of G that were ignored in estimation due to insufficient data
            Gregions = table1.region~=52 & table1.region~=88 & table1.region~=94; 
            rows = rows & Gregions;
        end
        
        T = table1(rows,:);

        % Get dummy vars (crop1Prev, t, region) from data table
        dummies = T{:,["crop1Prev" "t" "region"]};
        % Get ethanol index (actual) from data table;
        ethindex = T{:,["z2" "z2cf"]};
        % Get field acres from data table
        acres = T{:,"acres"};
        
        % Set up holders for iteration results
        holder(i).baseline = zeros(N,12);
        holder(i).cf_local = holder(i).baseline;
        holder(i).cf_full = holder(i).baseline;
        
        
        for ii=1:N
            % Draw random parameter values from distribution
            b_i = (mvnrnd(b,V))';
            %b_i = b;
            
            % Get time fixed effects parameters
            % delta = [delta_c delta_s]
            delta = [b_i(est.outcome==1 & strcmp(est.variable,"t")) b_i(est.outcome==5 & strcmp(est.variable,"t"))];
            % Call function to compute cf values of fixed effects,
            % delta_tilde
            delta_tilde = [fitRFS(delta(:,1)) fitRFS(delta(:,2))];
            
            % Create cf b_i vector
            b_i_cf = b_i;
            b_i_cf(est.outcome==1 & strcmp(est.variable,"t")) = delta_tilde(:,1);
            b_i_cf(est.outcome==5 & strcmp(est.variable,"t")) = delta_tilde(:,2);
            
            % Find areas of interest using functions below
            holder(i).baseline(ii,:) = expectedArea(b_i,est,ethindex(:,1),dummies,acres);
            holder(i).cf_local(ii,:) = expectedArea(b_i,est,ethindex(:,2),dummies,acres);
            holder(i).cf_full(ii,:) = expectedArea(b_i_cf,est,ethindex(:,2),dummies,acres);

            %--- Update waitbar ---%
            iter = (j-1)*length(sizegroups)*N + (k-1)*N + ii;
            waitbar(iter/(length(LRRs)*length(sizegroups)*N),f,"Done with LRR " + lrr + ", iteration " + ii);
            
        end
    end
end
toc
%% --- Add in sameCrop observations ---%
rows = table1.sameCrop==1 & (~strcmp(table1.lrrsym,"G") | (strcmp(table1.lrrsym,"G") & Gregions)) & (~strcmp(table1.lrrsym,"O") & ~strcmp(table1.lrrsym,"P") & ~strcmp(table1.lrrsym,"R")) & table1.acres<10;
T = table1(rows,:);
acres_c = T{:,"acres"}.*(T{:,"crop"}==1);
acres_s = T{:,"acres"}.*(T{:,"crop"}==5);
acres_o = T{:,"acres"}.*(T{:,"crop"}==9);

total_baseline = [sum(acres_c) sum(acres_s) sum(acres_o) sum(acres_c) 0 0 0 sum(acres_s) 0 0 0 sum(acres_o)];
total_cf_local = total_baseline;
total_cf_full = total_baseline;

for i=1:length(LRRs)
    total_baseline = total_baseline + holder(i).baseline;
    total_cf_local = total_cf_local + holder(i).cf_local;
    total_cf_full = total_cf_full + holder(i).cf_full;
end

% Convert to million acres/year
total_baseline = total_baseline/(11*10^6);
total_cf_local = total_cf_local/(11*10^6);
total_cf_full = total_cf_full/(11*10^6);

change_local = total_baseline - total_cf_local;
change_full = total_baseline - total_cf_full;

std_local = [std(total_baseline)' std(total_cf_local)' std(change_local)'];
std_full = [std(total_baseline)' std(total_cf_full)' std(change_full)'];


%% --- FUNCTIONS ---%
% Function fitRFS: return counterfactual timed fixed effect parameter
% values
function delta_tilde = fitRFS(delta)
RFS =    [0 0 0 0 0 4 4.7 9   10.5 12  12.6 13.2 13.8 14.4 15  15  15  15]';
RFS_cf = [0 0 0 0 0 4 4.7 5.4 6.1  6.8 7.4  7.5  7.5  7.5  7.5 7.5 7.5 7.5]';

X = [ones(18,1) RFS];
B = X\delta;

e = delta-X*B;
h = zeros(length(delta),1);
H = zeros(length(delta));
for i = 1:length(delta)
    h(i) = X(i,:)/(X'*X)*X(i,:)';
    H(i,i) = e(i)^2/(1-h(i))^2;
end

HC3 = (X'*X)\X'*H*X/(X'*X);

phi = B(2) + HC3(2,2)*trnd(16);

delta_tilde = delta-phi*(RFS-RFS_cf);
end

% Function croppatterns: return total Mha for crop totals and crop
% sequences
function [croptot, cropseq] = croppatterns(cropareas,totalarea)
Mha = cropareas*(0.404686/10^6)/11;
croptot = [Mha(:,1:2) totalarea-Mha(:,1)-Mha(:,2)];
cropseq = [Mha(:,3:4) Mha(:,1)-Mha(:,3)-Mha(:,4) Mha(:,5:6) Mha(:,2)-Mha(:,5)-Mha(:,6) Mha(:,1)-Mha(:,3)-Mha(:,5) Mha(:,2)-Mha(:,4)-Mha(:,6)];
cropseq(:,9) = croptot(:,3)-cropseq(:,7)-cropseq(:,8);
end

% Function feValues: look up (random) fixed effect values in b_i array,
% using estimates table
function fe = feValues(b_i,est,crop,varstring,dummyvec)
b_i_var = b_i(est.outcome==crop & strcmp(est.variable,varstring));
[~,loc] = ismember(dummyvec,est{strcmp(est.variable,varstring),"value"});
fe = b_i_var(loc);
end

% Function transitionProb: Compute transition probabilities given parameter
% values and data
function [p_xc, p_xs, p_xo] = transitionProb(cons,b_z,z,fe_crop1Prev,fe_t,fe_region)
V_corn = cons(1) + b_z(1)*z + fe_crop1Prev(1) + fe_t(:,1) + fe_region(:,1);
V_soy = cons(2) + b_z(2)*z + fe_crop1Prev(2) + fe_t(:,2) + fe_region(:,2);

p_xc = exp(V_corn)./(1+exp(V_corn)+exp(V_soy));
p_xs = exp(V_soy)./(1+exp(V_corn)+exp(V_soy));
p_xo = 1./(1+exp(V_corn)+exp(V_soy));
end

% Function totalArea: Compute total areas of interest given parameter
% values and data
function results = expectedArea(b_i,est,z,dummies,acres)
b_i_z = [b_i(est.outcome==1 & strcmp(est.variable,"z2")) b_i(est.outcome==5 & strcmp(est.variable,"z2"))];
b_i_cons = [b_i(est.outcome==1 & strcmp(est.variable,"_cons")) b_i(est.outcome==5 & strcmp(est.variable,"_cons"))];

fe_cornPrev = [feValues(b_i,est,1,"crop1Prev",1) feValues(b_i,est,5,"crop1Prev",1)];
fe_soyPrev = [feValues(b_i,est,1,"crop1Prev",5) feValues(b_i,est,5,"crop1Prev",5)];
fe_otherPrev = [feValues(b_i,est,1,"crop1Prev",9) feValues(b_i,est,5,"crop1Prev",9)];
fe_t = [feValues(b_i,est,1,"t",dummies(:,2)) feValues(b_i,est,5,"t",dummies(:,2))];
fe_region = [feValues(b_i,est,1,"region",dummies(:,3)) feValues(b_i,est,5,"region",dummies(:,3))];

[p_cc,p_cs,p_co] = transitionProb(b_i_cons,b_i_z,z,fe_cornPrev,fe_t,fe_region);
[p_sc,p_ss,p_so] = transitionProb(b_i_cons,b_i_z,z,fe_soyPrev,fe_t,fe_region);
[p_oc,p_os,p_oo] = transitionProb(b_i_cons,b_i_z,z,fe_otherPrev,fe_t,fe_region);

P_c = (p_os.*(p_sc-p_oc) + p_oc.*(1-p_ss+p_os))./((1-p_ss+p_os).*(1-p_cc+p_oc) - (p_sc-p_oc).*(p_cs-p_os));
P_s = (p_oc.*(p_cs-p_os) + p_os.*(1-p_cc+p_oc))./((1-p_ss+p_os).*(1-p_cc+p_oc) - (p_sc-p_oc).*(p_cs-p_os));
P_o = 1-P_c-P_s;

acres_c = acres.*P_c;
acres_s = acres.*P_s;
acres_o = acres.*P_o;

acres_cc = acres_c.*p_cc;
acres_cs = acres_c.*p_cs;
acres_co = acres_c.*p_co;
acres_sc = acres_s.*p_sc;
acres_ss = acres_s.*p_ss;
acres_so = acres_s.*p_so;
acres_oc = acres_o.*p_oc;
acres_os = acres_o.*p_os;
acres_oo = acres_o.*p_oo;

m = [acres_c acres_s acres_o acres_cc acres_cs acres_co acres_sc acres_ss acres_so acres_oc acres_os acres_oo];
results = sum(m);
end
